
#include "expo_fusion.h"

expo_fusion::expo_fusion()
{
}

int expo_fusion::process(Mat &src)
{
	clock_t t0 = clock();

	vector<Mat> w(NIMAGE);	//weight map
	vector<Mat> g(NIMAGE);	//gamma corrected image

	//------------------------------------------------------------------------------
	//	Desc: synthesis several gamma corrected image
	//------------------------------------------------------------------------------
	gamma_correc(src, g.at(0), 2.2f);		//gamma 2.2
	gamma_correc(src, g.at(1), 1.0f/2.2f);	//gamma	0.45 (1/2.2=0.45)
	gamma_correc(src, g.at(2), 1.0f);		//gamma 1 (original)

	clock_t t1 = clock();
	printf("gamma correct %d\n",t1-t0);

	//------------------------------------------------------------------------------
	//	Desc: synthesis and normalize weight maps
	//------------------------------------------------------------------------------
	for(int i=0; i<NIMAGE; i++)
		getWeightMap(g.at(i), w.at(i));

	normalize_weight(w);

	clock_t t2 = clock();
	printf("synthesis weight maps %d\n",t2-t1);

	//------------------------------------------------------------------------------
	//	Desc: build laplacian pyramid for each image
	//------------------------------------------------------------------------------
	vector<vector<Mat>> pyramid_laplacian(NIMAGE);

	for(int i=0; i<NIMAGE; i++)
	{
		pyramid_laplacian.at(i).resize(NPYRAMID);
		build_laplacian_pyramid(g.at(i), pyramid_laplacian.at(i));
	}

	clock_t t3 = clock();
	printf("build laplacian pyramid %d\n",t3-t2);

	//------------------------------------------------------------------------------
	//	Desc: build gaussian weight pyramid each weight map
	//------------------------------------------------------------------------------
	vector<vector<Mat>> pyramid_gaussian_weight(NIMAGE);
	for(int i=0; i<NIMAGE; i++)
	{
		pyramid_gaussian_weight.at(i).resize(NPYRAMID);
		buildPyramid(w.at(i), pyramid_gaussian_weight.at(i), NPYRAMID-1);
	}

	clock_t t4 = clock();
	printf("build gaussian weight pyramid %d\n",t4-t3);

	//------------------------------------------------------------------------------
	//	Desc: pyramid blending
	//------------------------------------------------------------------------------
	for(int i=0; i<NPYRAMID; i++)
	{
		for(int y=0; y<pyramid_laplacian.at(0).at(i).rows; y++)
		{
			for(int x=0; x<pyramid_laplacian.at(0).at(i).cols; x++)
			{
				for(int k=0; k<NIMAGE; k++)
				{
					pyramid_laplacian.at(k).at(i).at<float>(y,x*3)	 *= pyramid_gaussian_weight.at(k).at(i).at<float>(y,x);
					pyramid_laplacian.at(k).at(i).at<float>(y,x*3+1) *= pyramid_gaussian_weight.at(k).at(i).at<float>(y,x);
					pyramid_laplacian.at(k).at(i).at<float>(y,x*3+2) *= pyramid_gaussian_weight.at(k).at(i).at<float>(y,x);
				}

				for(int k=1; k<NIMAGE; k++)
				{
					pyramid_laplacian.at(0).at(i).at<float>(y,x*3)   += pyramid_laplacian.at(k).at(i).at<float>(y,x*3);
					pyramid_laplacian.at(0).at(i).at<float>(y,x*3+1) += pyramid_laplacian.at(k).at(i).at<float>(y,x*3+1);
					pyramid_laplacian.at(0).at(i).at<float>(y,x*3+2) += pyramid_laplacian.at(k).at(i).at<float>(y,x*3+2);
				}
			}
		}
	}

	clock_t t5 = clock();
	printf("pyramid blending %d\n",t5-t4);
	
	//------------------------------------------------------------------------------
	//	Desc: reconstruction
	//------------------------------------------------------------------------------
	Mat recon;
	reconstruction(pyramid_laplacian.at(0), recon);

	clock_t t6 = clock();
	printf("reconstruction %d\n",t6-t5);

	//------------------------------------------------------------------------------

	return 0;
}

//------------------------------------------------------------------------------------------
//	Desc: use Gamma correction to synthesis several input images
//------------------------------------------------------------------------------------------
void expo_fusion::gamma_correc(Mat &src, Mat&dst, float gamma)
{
	dst.create(src.size(), CV_32FC3);

	if(gamma == 1.0f)
	{
		//------------------------------------------------------------------------------------------
		//	Desc: using pow() is time consuming, so avoid call pow() as possible
		//------------------------------------------------------------------------------------------
		for(int y=0; y<src.rows; y++)
		{
			for(int x=0; x<src.cols; x++)
			{
				dst.at<float>(y, x*3)   = src.at<unsigned char>(y, x*3  )/255.0f;
				dst.at<float>(y, x*3+1) = src.at<unsigned char>(y, x*3+1)/255.0f;
				dst.at<float>(y, x*3+2) = src.at<unsigned char>(y, x*3+2)/255.0f;
			}
		}
	}
	else
	{
		for(int y=0; y<src.rows; y++)
		{
			for(int x=0; x<src.cols; x++)
			{
				dst.at<float>(y, x*3)   = pow(src.at<unsigned char>(y, x*3  )/255.0f, gamma);
				dst.at<float>(y, x*3+1) = pow(src.at<unsigned char>(y, x*3+1)/255.0f, gamma);
				dst.at<float>(y, x*3+2) = pow(src.at<unsigned char>(y, x*3+2)/255.0f, gamma);
			}
		}
	}
}

//------------------------------------------------------------------------------
//	Desc: Generate a weight map based on several quality measures.
//------------------------------------------------------------------------------
int expo_fusion::getWeightMap(Mat &src, Mat &weight)
{
	clock_t t0 = clock();

	Mat gray;
	cvtColor(src, gray, CV_BGR2GRAY);
	
	clock_t t1 = clock();
	printf("	CV_BGR2GRAY %d\n",t1-t0);

	//------------------------------------------------------------------------------
	//	Desc: well-exposedness measurement
	//------------------------------------------------------------------------------
	Mat expo(src.size(), CV_32FC1);
	for(int y=0; y<src.rows; y++)
	{
		for(int x=0; x<src.cols; x++)
		{
			float luma = gray.at<float>(y, x);
			float sigma = 0.2f;

			expo.at<float>(y, x) = exp(-0.5f * (float)pow(((luma-0.5f)/sigma), 2));
		}
	}

	clock_t t2 = clock();
	printf("	exposure mesure %d\n",t2-t1);

	//------------------------------------------------------------------------------
	//	Desc: contrast measurement
	//------------------------------------------------------------------------------
	Mat contr = Mat::zeros(src.size(), CV_32FC1);
	for(int y=1; y<src.rows-1; y++)
	{
		for(int x=1; x<src.cols-1; x++)
		{
			 float contrast =  -gray.at<float>(y, x)*4.0f + 
								gray.at<float>(y, (x+1)) + gray.at<float>(y, (x-1)) +
								gray.at<float>((y+1), x) + gray.at<float>((y-1), x);

			 contr.at<float>(y, x) = abs(contrast);
		}
	}

	clock_t t3 = clock();
	printf("	contrast mesure %d\n",t3-t2);

	//------------------------------------------------------------------------------
	//	Desc: saturation measurement
	//------------------------------------------------------------------------------
	Mat satur = Mat::zeros(src.size(), CV_32FC1);
	for(int y=0; y<src.rows; y++)
	{
		for(int x=0; x<src.cols; x++)
		{
			float b = src.at<float>(y,x*3);
			float g = src.at<float>(y,x*3+1);
			float r = src.at<float>(y,x*3+2);

			float maximum = max(max(b, g), r);
			float minimum = min(min(b, g), r);

			if(maximum == minimum)
				satur.at<float>(y,x)=0;
			else
			{
				float sum = maximum + minimum;
				float diff = maximum - minimum;

				if(sum<=1.0f)
					satur.at<float>(y,x) = diff/sum;
				else
					satur.at<float>(y,x) = diff/(2.0f-sum);
			}
		}
	}

	clock_t t4 = clock();
	printf("	saturation mesure %d\n",t4-t3);

	//------------------------------------------------------------------------------
	//	Desc: merge
	//------------------------------------------------------------------------------
	weight.create(src.size(), CV_32FC1);

	float satur_weight = 1.0f;
	float contr_weight = 1.0f;
	float expo_weight = 1.0f;

	if((satur_weight == 1.0f) && (contr_weight == 1.0f) && (expo_weight == 1.0f))
	{
		//------------------------------------------------------------------------------------------
		//	Desc: using pow() is time consuming, so avoid call pow() as possible
		//------------------------------------------------------------------------------------------

		for(int y=0; y<src.rows; y++)
			for(int x=0; x<src.cols; x++)
				weight.at<float>(y, x) = satur.at<float>(y, x) * contr.at<float>(y, x) * expo.at<float>(y, x);		
	}
	else
	{
		for(int y=0; y<src.rows; y++)
		{
			for(int x=0; x<src.cols; x++)
			{
				weight.at<float>(y, x) = pow(satur.at<float>(y, x), satur_weight) * 
										 pow(contr.at<float>(y, x), contr_weight) *
										 pow(expo.at<float>(y, x),  expo_weight);
			}
		}
	}

	clock_t t5 = clock();
	printf("	weight merge %d\n",t5-t4);
	printf("\n");

	return 0;
}

//------------------------------------------------------------------------------
//	Desc: normalize all weight maps such that weight values sum to one
//		  at each pixel.
//------------------------------------------------------------------------------
void expo_fusion::normalize_weight(vector<Mat> &m)
{
	int nMat = (int)m.size();

	vector<float> r(nMat);

	for(int y=0; y<m.at(0).rows; y++)
	{
		for(int x=0; x<m.at(0).cols; x++)
		{
			bool isAllZero = true;

			for(int i=0; i<nMat; i++)
			{
				r.at(i) = m.at(i).at<float>(y,x);

				if(r.at(i) != 0)
					isAllZero = false;				
			}

			if(isAllZero)
			{
				//------------------------------------------------------------------------------
				//	Desc: if all weights are zero, normalized weights average value.
				//		  (but the correctness haven't been confirm)
				//------------------------------------------------------------------------------

				for(int i=0; i<nMat; i++)
					m.at(i).at<float>(y,x) = 1.0f/nMat;
			}
			else
			{
				float sum = 0.0f;

				for(int i=0; i<nMat; i++)
					sum += r.at(i);

				for(int i=0; i<nMat; i++)
					m.at(i).at<float>(y,x) = r.at(i) / sum;
			}
		}
	}
}

//------------------------------------------------------------------------------------------
//	Desc: Lossless Image Decomposition, please refer to "The Laplacian Pyramid as a Compact 
//		  Image Code Burt, P., and Adelson, E. H., IEEE Transactions on Communication, 
//		  COM-31:532-540 (1983)."
//------------------------------------------------------------------------------------------
void expo_fusion::build_laplacian_pyramid(Mat &src, vector<Mat> &pyramid_laplacian)
{
	vector<Mat> pyramid_gaussian;
	buildPyramid(src, pyramid_gaussian, NPYRAMID-1);

	for(int i=0; i<NPYRAMID-1; i++)
	{
		Mat pUp;
		pyrUp(pyramid_gaussian.at(i+1), pUp);

		pyramid_laplacian.at(i) = pyramid_gaussian.at(i)- pUp;
	}

	pyramid_laplacian.back() = pyramid_gaussian.back().clone();

	
}

//------------------------------------------------------------------------------------------
//	Desc: Lossless Image reconstruction, please refer to "The Laplacian Pyramid as a Compact 
//		  Image Code Burt, P., and Adelson, E. H., IEEE Transactions on Communication, 
//		  COM-31:532-540 (1983)."
//------------------------------------------------------------------------------------------
void expo_fusion::reconstruction(vector<Mat> &pyramid_laplacian, Mat &recon)
{
	pyrUp(pyramid_laplacian.at(NPYRAMID-1), recon);

	for(int i=NPYRAMID-2; i>=1; i--)
	{
		recon += pyramid_laplacian.at(i);
		pyrUp(recon, recon);
	}

	recon += pyramid_laplacian.at(0);
	imshow("recon", recon);

	//imwrite("C:/Users/kuenbutsai/Desktop/New folder/result.bmp", recon*255.0f);
}